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(g) Vorrichtung fur SPECT-Untersuchungen 

(g) Die Erfindung betrifft eine Vorrichtung sowie ein Ver- 
fahren zur Tomograph ie, insbesondere zur Einzelgamma- 
quanten-Tomographie (SPECT). 

Die Vorrichtung zur Durchfuhrung eines tomographi- 
schen Verfahrens, insbesondere zur Durchfuhrung einer 
Einzelphotonen-Tomographie, weist einen Multi-Loch- 
Kollimator und einen Detektor zur Erfassung von Gamma- 
quanten oder Photonen auf, die durch den Multi-Loch- 
Kollimator hindurchtreten. 

Beim Verfahren zur Durchfuhrung des tomographischen 
Verfahrens ist der Abstand zwischen dem Objekt und dem 
Multi-Loch-Kollimator kleiner gewahit als der Abstand 
zwischen dem Multi-Loch-Kollimator und der Oberflache 
des Detektors. 

Es kann mit hoher Ortsauflosung und Sensitivitat das ge- 
■ wunschte Ergebnis ermittelt werden. 
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Beschreibung 

[0001] Die Erfindung betrifft eine Vbrrichtung sowie ein Verfehren zur Tomogr^hie. insbcsondere zur Einzelgamma- 

s fS" SriSphoS'o^etTomogra bezieht sich auf ein Verf ahren nebst zugehorigen Votrichtungen zur dreidi- 

" mensionalen Darstellungen von Radiopharmaka, die in ein Objekt gebracht wurden. Als Objekt konnen Menscjen oder 
Tiere vorgesehen sein. Die in das Objekt gebrachten Radiopharmaka senden Gamniaquanten. Die Oaininaquanten wer-, 
den von der Vorrichtung erfasst und ausgewertet. Als Ergebnis der Auswenung wird die Lage, also die raunUicte >ferter- 
lung der Radiopharmaka im Objekt erhalten. Die Lage der Radiopharmaka erlaubt wiedemm Ruckschlusse aut das Ob- 

10 iekt so zumBeispiel auf eine Verteilung von Gewebeim Objekt. 

[0003] Eine betanntc Vorrichtung zur Duichfiihrung einer Einzelgammaquanten-Tomographie umtass eine Gamina- 
kamera undeinen vorgeschalteten KoUimator. Beim KoUimatorhandelt es sich im aUgenieinen um eine Bleiplatte nut ei- 
ner Vielzah) von senl^ht durch die Platte fUhrenden Kanalen. Durch das \brsehen der Kanale wird sichergestellt, dass 
einerseits nur senkrecht einfallende Gammaquanten erfasst werden und andererseits eme ortsauflosende Messung mog- 

15 lich ist Die Kamera wird samt KoUimator um das Objekt herum verfahren. Hierdurch erhalt man eine Vielzahl von Orts- 
informationen. Es handelt sich dabei um sogenannte Projektionsaufiiahmen. Aus den gewonnenen Ortsinformationen 
rundum das Obiekt herum kannanschlieBend die Lage der Radiopharmaka im Objekt ermittelt werden. 
[00041 Um von den Gammaquanten verursachte Su«ustrahlung ausblenden zu konnen, wird im aUgememen n«:h eine 
Energieinformation benotigt. Die Kamera ist daher in der Regel so beschaffen, dass mit ihr sogleich die Energie der em- 

•>0 fallenden Gammaquanten ermittelt werden kann. ,. . ui Co 

[0005] Streustrahlung weist grundsatzlich eine geringere Energie auf im Vergleich zur eigenthchen Messsu^ahlung Es 
kann so Streustrahlung ausgeblendet werden, indem Gammaquanten mit niedriger Eneigie nicht berucksichugt werderi. 
Eine Obeigrenze der Energie der Gammaquanten festzulegen kann ebenfalls von Interesse sein, um Hmteigrundsttah- 

luDK ausblenden zu konnen. • t- u ■ „ ^« 

[0006] Das vorbeschriebene Verfahren bzw. die vorbeschriebene Vorrichtung gehort zum aUgemeinen Fachwissen, da 
solche Verfahren und \foirichlungenseitmehr als dreiBigJahreneingeselzt werden. ^^^r. , „ t 

[0007] Die Einzclgammaquanlcn-Tomographic (SPECT) und die Positroncn-Emissionstomographic (PET) stcUcn In- 
strumente zur quanUtativen Darstellung raumlicher Radiotracer-Verteilungen in vivo dar. AuBer in der Humanniedizin 
lassen sich diese Verfahren in der pharmakologischen und praklinischen Forschung zur Entwicklung und Evaluation 
neuartiger Tracer- Verbindungen einseizen. Wahrend in der PET heute diverse Systeme zur Untersuchung klemer Labor- 
tiere zur Verfugung stehen, hat es im Beieich der SPECT entsprechende Entwicklungen bisher nicht oder in nur unzuret- 
chendem MaBe gegeben, und das obwohl TC-99m und M23-markierte Radiopharmaka m der Nuklearmedizin eme un- 
glach hehere Bedeutung haben ak die PET-Nuklide. ... , . . u tt u a 

[0008] Mit einem hochauflosenden und iKxrhsensitiven Tier-SPECT ergabe sich fur die praklinische Forschung der 
Vorteil eines tierschonenden Verfahrens, mit dem sich aussagekrafuge Studien dynamisch und wiederholbar an emem In- 
dividuum durchfuhren lieBen. Dies wird dadurch begunstigt, dass sich bei den oben genannten Radioisot^n extrem 
hohe spezifische Aktivitaten erzielen lassen (ca. Faktor 100 gegenuber PET-Nuklideri) welche tiir storungsfeie^n vivo 
Messungen unerlassUch sind (geringe Massendosis). Hierzu sind begleitende Entwicklungen entsprechender Markie- 

[OOoS^^S *e O^SuflSung gegenuber dem eingangs genannten Stand der Technik zu verbessem wird ein l^ch- 
koUimator bei der Einzelgammaquanten-Tomographie eingesetzt. Ein lxx:h-Komiiiator zeichnet ^'^^^"^^^^^"^^^ 
nes Loch aus. durch das die Gammaquanten hindurchtreten. Befindet sich das Objekt naher am I^h-KoUimator als die 
Oberflache einer Gammakamera bzw. eines Detektors, so wird hierdurch schlieBlich eine vergroBerte Ortsauflosung er- 
reicht Durch den Loch-KoUimator hindurch treten die Gammaquanten nicht ausschbeBhch senkrecht. Stattdessen tr«ea 
45 sie kegelfdrmig ein und wieder aus. Da der hinter dem lx,ch-Kollimator liegende Kegel groBer ist als der Kegel vor dem 
Loch-Kollimator, wird im Ergebnis eine Veibesserung der Ortsauflosung im Vergleich zum eingangs genannten Stand 

fooSf 'T^^^m Loch-KoUimator ist eine kleine Durchtriltsoffhung bzw. ein kleines Loch vorzusehen, durch das die 
Gammaquanten hindurchtreten. um eine gute Ortsauflosung zu erhalten. Je kleiner ein Loch j^^oc ^st, desto we^^^^^^ 
50 Gammaquanten treten durch dieses Loch hinduich. Mit kleiner weidendem l^h sinkt daher nachieilhaft die Sensmyitat 
der Vorrichmng. SenativitSt ist definiert als das Verhatnis von gemessener Zahkate zur im Objekt vothandenen Aktivi- 

[Mll] Wird die Sensitivitat zu gering, so ist schlieBlich die Durchfiihrung raner Knzelgammaquanten-Tomographie 

55 [oSir AiS'Se der Erfindung ist die Schaffung einer Vorrichtung nebst zugehorigem Verfahren der eingangs genann- 
ten Art, mit der hochauflosend und hochsensitiv gemessen werden kann. . . u 
[00131 Die Aufgabe der Erfindung wird durch eine Vorrichtung mit den Merkmalen des ersten Anspruchs sowie durch 
ein Verfahren mit den Merkmalen des Nebenanspruchs gelost. Vorteilhafte Ausgestaltungen ergeben sich aus den Unter- 

60 Sff*Die anspruchsgemaBe VorrichUmg umfasst einen Muhi-7x.ch-Kollimator nebst einem Detektor zur Erfa.ssung 
der Gammaquanten. die durch den Multi-Loch-KolHmator hindurchtreten. Der KoUimator wejst also eine Mehrz^l an 
Durchtrittsoffnungen auf. In einer Ausgestaltung der Erfindung ist der Detektor so beschaften, dass dieser auch die Ener- 
gie der aufueffenden Gammaquanten zu bestimmen vermag. r.^,T3;„ 
100151 Da der KoUimator mehrere Locher aufweist, steigt die Sensitivitat der VDmchtung entsprechend an. Der Ein- 

65 satz eines Loch-KoUimaiors hat gegenUber dem Einsatz solcher Kollimatoren, mil denen nur senkrecht einfallende 
Strahlen erfasst werden, den Vorteil der hohen Ortsauflasung. Somit steht eine Vorrichtung mit guier Ortsauflosung und 
guter Sensitivitat zur Verfugung. „ . . , , • t u ly ir ,™ i?« 

[0016] Wahrend des Betriebes der Vorrichtung befindet sich das Objekt naher am Multi-Loch-Kolhmator als die Ka- 



25 



30 



35 



40 



SIXXID: <DE. 



.10142421A1J_> 



101 42 421 A 1 



mera- bzw. Detektoroberflache, uni eine gute Ortsauflosung zu erzielen, Bei der \brrichtung befindet sich die Halterung 
fur das Objekt (Patientenliege) daher naher am Multi-Lx)ch-Kollimalor als die Kamera bzw. der Detektor. 
[0017] Die Abslande der einzelnen Durchtritle oder Locher im Multi-Loch-Kollimator sind vorzugsweise so gewahlt, 
dass sich die Kegel, die auf der Kamera auftreften, allenfalls teilweise iiberschneiden. Zur Erzielung einer guten Ortsauf- 
losung sowie guter Sensitivitat ist es vorteilhaft, Uberschneidungsbereiche zuzulassen. Diese betragen in einer Ausge- 5 
staltung der Erfindung nicht mehr als 30%, vorzugsweise bis zu 70% der Gesamtfiache eines Kegels, der durch die Gam- 
nfiaquanten gebildet wird, die durch ein Loch des Multi-Loch-Kollimators hindurchireten. 

[d018] In der herkommlichen Loch-Tomographie (Pinhole-Tomographic) liegt das Lochzentruni auf der Miltensenk- 
rechten des Deteklors. Femer sleht die Loch-Achse, d. h. die Symmetrieachse des KoUimatorlochs, senkrecht zum De- 
tektor. Es wird dann ein Rekonstruktionsverfahren verwendet, bei dem davon ausgegangen wird, dass das Gamiiiaquant, lo 
welches voni Zentrum des Objekts auf die Kamera fallL, mit der Kameraoberflache einen rechten Winkel bildet. Die 
Grundflache des Kegels, die auf der Kamera abgebildet wird, ist dann grundsatzlich kreisformig, 

[0019] Diese Situation ist. beini Hinsaty. eines Multi-T^ch-Kollimators vielfach nicht der Fall. F^s wird daher in einer 
Ausgestaltung der Erfindung ein Rekonstruktionsverfahren bereit gestellt, welches die abweichenderi Bedingungen be- 
riicksichtigt. Fallt ein vom Zentrum des Objekts stammendes Gammaquant nicht mehr senkrecht auf die Oberflache der 15 
Kamera bzw. des Deteklors auf, so wird kein kreisfdrmiger Kegel (idealisierter Zustand) auf der Kameraoberflache ab- 
gebildet. Stattdessen bildet sich der Kegel grundsatzlich in Form einer Ellipse auf der Kamem ab. ErfindungsgemaB wird 
dieses Problem gel5st, indem ein iterativcs Rekonstruktionsverfahren eingesetzt wird. Ausgangspunkt des iterativen Re- 
konstruktionsverfahrens ist eine angenommene Verteilung im Objekt und zwar in der Regel eine raumliche Verteilung. 
Es wird dann berechnet, welches Messergebnis die angenommene Verteilung erzielen wiirde. Das berechnete Ergebnis 20 
wird mit dem tatsachlich gemessenen verglichen. AnschlieBend wird eine neue geanderte Verteilung genommen. Wie- 
derum wird das auf der Kamera abgebildete Ergebnis dieser neuen Verteilung berechnet. Es wird emeut verglichen. Es 
wird festgestellt, ob die neue Verteilung dem gemessenen Ergebnis besser entspricht. Auf diese Weise wird nach Durch- 
fuhrung hinreichend vieler Schritte eine Verteilung ermittelt, deren berechnetes Ergebnis mit dem tatsachhchen Ergebnis 
(Messergebnis) hinreichend gut ubereinstimmt. Das iterative Rekonstruktionsverfahren ist insbesondere dann beendet, 25 
wenn mil einer vorgegebenen Gehauigkeil das berechnete Ergebnis mil dem gemessenen ubereinsliimiit. Das Iterations- 
vcrfahrcn umfasst also cine sogcnanntc Vorwartsprojcktion, also die Bcrcchnung des Ergcbnisscs cincr angcnommcncn 
Verteilung. 

[0020] Das iterative Verf ahren weist femer den Vorteil auf, dass Uberlappungsbereiche der auf der Kamera- bzw. De- 
tektoroberflache abgebildeten Kegel berechnet und mit dem tatsachlichen Eigebnis verglichen werden konnen. Es ist da- 30 
her auch aus diesem Grund anderen Rekonstruktionsverfahren vorzuziehen. Es ist also so moglich, Uberlappungsberei- 
che zuzulassen und so zu guten Ortsauflosungen zu gelangen. 

[0021] In einer weitercn Ausgestaltung der Erfindung umfasst der Multi-Loch-KoUimaior eine Platte, die aus Wolfram 
und Iridium gefertigt ist. Diese Materialien besitzen einen besseren Schwachungskoeffizienten gegenuber Gammaquan- 
ten im Vergleich zu Blei, Iridium ist unter den genannten Materialien am besten geeignet, um Geunmaquanten zu schwa- 35 
chen. Iridium ist jedoch sehr teuer. Daher wird aus Kostengriinden Wolfram an den Stellen eingesetzt, bei denen die An- 
fordemngen an das Schw^achungsverhalten weniger groB sind. Aus Iridium werden die Teile der Platte gefertigt, bei de- 
nen die Anforderung^ an die Schwachung von Gammaquanten besonders groB sind. Es handelt sich dabei insbesondere 
um die Bereiche der Platte, die an die Locher angrenzen. 

[0022] Ein toch in der Platte mundet vorteilhaft von beiden Seiten trichterformig in die Platte ein. Hier sind die An- 40 
forderungen an die Schwachung besonders groB und zwar insbesondere an der Lochwandung. Daher sind die Trichter- 
wande bevorzugt aus Iridium gefertigt- Die Platte ist dann typischerweise 3 bis 10 mm dick. 

[0023] Gammaquanten, die vom Inneren des Objektes ausgehen, werden in der Regel gewebeabhangig geschwacht, 
GemaB Stand der Technik wird zur Beriicksichtigung dieser Schwachung bei der Auswertung von einem homogenen 
Schwachungskoeffizienten ausgegangen, der dem Schwachungskoeffizienten von Wasser entspricht. Daneben ist die 45 
Schwachung von den Umrissen des Objektes abhangig. In einer Ausgestaltung der Erfindung wird im Rahmen der Aus- 
wertung der auBere Umriss des Objektes ermittelt und die Schwachung in Abhangigkeit vom Umriss kalkuliert. Auf 
diese Weise werden weiter verbesserte Ergebnisse erhalten. 

[0024] MaB fur den auBeren Umriss des Objektes ist die Compton-Streustrahlung. In einer Ausgestaltung der Erfin- 
dung wird daher die Compton-Streustrahlung zum Beispiel in einem sogenannten Compton-Fenster gemessen. Im Re- 50 
konstruktionsverfahren wird die Compton-Streustrahlung beriicksichtigt und hieraus der Umriss des Objektes. ermittelt. 
[0025] Trifft ein Gammaquant auf eine Kamera bzw. den Detektor auf, so wird der Ort des Auftreffens mit einer ka- 
mera- Oder detektortypischen Ungenauigkeit gemessen. In einer weiteren Ausgestaltung des Verfahrens wird bei der Vor- 
wartsprojcktion, die dem iterativen Rekonstruktionsverfahren zugrunde liegt, die Abbildungseigenschaft, also die ka- 
mera- oder detektortypische Ungenauigkeit bei der Auswertung beriicksichtigt. Wiederum gelingt die Beriicksichtigung 55 
der Messungenauigkeit durch ein Iterationsverfahren der vorgenannten Art zuverlassig. 

[0026] Befindet sich eine im Objekt befindliche strahiende Quelle verhaltnismaBig weit von dem Multi-Loch-Kollima- 
tor entferat (also in einem Bereich des Objekts, der besonders weit weg vom Kollimator liegt), so ninmit die Sensitivitat 
ab. In einer Ausgestaltung des Rekonstruktionsverfahrens wird bei der Vorwartsprojcktion diese abnehmende Sensitivi- 
tat herucksichtigt. 60 
[0027] Das Abbildungsverhalten der Kamera bzw. des Deteklors hangt ebenfaUs von der Entfernung ab, die zwischen 
der slrahlenden Quelle und dem Multi-Ix)ch-KoIlimator besteht. Dieses sich entfernungsabhangig andemde Abbildungs- 
verhalten wird in einer Ausgestaltung des Verfahrens ebenfalls iterativ berucksichtigt. 

[0028] Nachfolgend werden einschlagige Programmteile flir ein Iterationsverfahren angegeben, welches die vorge- 
nannten erfindungsgemaBen Schritte abzuarbeiten vermag. Die Programme umfassen die nachfolgend genannten Ein- 65 
gabe-Parameter. Ferner werden typische Werte solcher Eingabeparameter angegeben. Der Begriff "Pinhole" wird syn- 
onym fiir den Begriff "Loch" (des Multi-Loch-Kollimators) verwendet. 

[0029] Soweit die Berechnung der Kontur bzw. des TJmrisses des Objekts betroflfen ist, wird die Berechnung der Ob- 
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jektkontur in einer Ausfiihrungsform in der "nuliten" Iteration vorgenommen. Eine eigentlich die Bildqualitat degradie- 
rende Eigenschaft der Comptonstreuung wird hierfiir ausgenutzt: die richtungsfalsch deiektierten Gamiftaquanien. 
[0030] Richtungsfalsch detektierte Gainmaquanten stellen einen Untergrund in den Projektionen dar, der nachhallig 
die Bildqualitat verschlechtert. Sie konnen jedoch auch von Nutzen sein, Sie bewirken, dass in praktisch alien kiinischen 
5 Fallen die gesamte Ausdehnung des Patienten in den Projektionen erscheint. Auch wenn sich ein Tracer, also ein Radio- 
pharmakon sehr spezifisch an einem eng umgrenzten Organ anlagem soUle, so scheinen dennoch auch Quanten vpn alien 
anderen zur Abbildung gebrachten Bereichen des Patienten zu stammen. Tatsachlich sind dies Quanten, die ihren Ur- 
sprung in dem eng umgrenzten Organ haben, aufgrund der Comptonstreuung jedoch scheinen sie den ganzen Patienten 
zu "illuiiiinieren", Dieser Urns land wird ausgenutzt, um die Objektkontur zu berechnen. 
10 [0031] Die Berechnung erfolgt mehrstufig: 

1. Erstellung "binarer" Projektionen. Sie stellen eine Vereinfachung der eigentlichen Projektionen dar, insofem, 
dass bei ihnen jeder Pixelinhalt groBer Null auf Kins gesetzt wird. 

Wesentlich datlur ist eine benutzergefuhrte Schwellsetzung, die nach Schliissigkeit die Projektion des eigentlichen 
15 Untarsuchungsobjektes, des Patienten also, vom Hintergrund separiert Die Berechnung der Schwelle orientiert sich 

an einem gemittelten Maximum, das aus samtlichen Projektionen gebildet wird. 

2. Riickprojektion der binaren Projektionen in den Objektraum. 

3. Anhand der "N^elfachheit" der Voxel (kleines meist kubisches Volumenelement), d. i. die Haufigkeit, mit der ein 
Voxel unter der jeweiligen KoUimatorgeometrie iiber alle Winkel von den Projektionen gesehen wird, setzt eine (fiir 

20 die jeweilige Geometric heuristisch emiittelte) Schwelle fest, welche Voxel in allererster Naherung zum Korperbin- 

nenraum gehoren. (Die Begrenzung des Korperbinnenraumes ist die Korperkontur.) 

4. Mehrmalige Faltung niit 3d-Faltungskem 

5. Wiederholung Punkt 3 

6. Zweimaliges Diirchlaufen von: 
25 a) 3d-Faltung 

b) Hinzunahiiie der Voxel, in die elwas hineingefallel wurde, zuiu Korperinnenrauiii. 



Eingabe-Aufforderung 

Projektionen (float): maus.prj 
Breite der Proj. [pix]: 266 
Schichten in Proj. [pix]: 193 
Anzahl der Winkel: 60 
Pixelgroesse in Proj. [mm]: 2 
Projektionen glaetten (nein = 0/ja = 1): 1 
Abstand RotAchse/Bildebene [mm]: 155 
Name des Aperturmusters (text): apertur.txt 

Dicke des Krislalls [mm]: 10 
Absorption des KristaUs [1/cm]: 1.173 
Intrinsische Aulloesung [mm]: 3.3 
Abs. Intensitaet des Gammas [%]: 90 



Bedeutung 

Die gemessenen Projektionen 
Transversale Dimension 
Axiale Dimension 
Anzahl der Projektionen 
GroBe der Projektionspixel 
Glattung der gemessenen Projektionen 
Abstand Bildebene zu Rotationsachse 
Name der Aperturdatei (siehe unten) 

Dicke des Szintillators 
Absorptionskoeffizient Kristall 
Intrinsische Autlosung der Kamera 
Abs. Intensitat der Ganunalinie 



Startverteilung (float): maus.sta 
Breite der Rekon. [pix]: 112 
Schichten der Rekon, [pix]: 262 
Voxelgroesse in Rekon. [mm]: 0.4 
Zahl der Iteration en: 30 
5^ Voiumen glaetten (nein = 0/ja = 1): 1 

Abbildungsfkt. voxelweise (nein = 0/ja = 1)?: 0 



Startverteilung der Rekonstruktion 
Transversale Dimension 
Axiale Dimension 
GroBe der Objektvoxel 
Anzahl der Iterationen 
Glattung der Volumendaten 

PSF voxelweise beriicksichtigen oder mit mittlerer PSF falten? 



Abw. vom Original (nein = 0/ia = 1)?; 0 Abweichung zum Original bestimmen? 

Stamnmame des Outputs: maus Namensstamm der Output-Dateien 

Speichem alle?: 5 Zwischenergebnisse speichem? 



Aufbau der Aperturdatei 

60 [0032] Fine Aperturdatei hat typischerweise folgenden Aufbau: 
7 
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[0033] Beim ersten Eintr^ handelt es sich um die Anzahl der Locher. Jede folgende Zeile beschreibt ein Pinhole, Die 
Bedeutung der Eintrage ist folgende: 



Spalte 


Bedeutung 




1 


x-Koordinate de-s Pinholemittelpunkts, d. h. Abstand des Pinholes von der Rotationsachse in mm 


5 


o 


^ y-Koord. des Pinholes, d. h. transversale Auslenkung des PH 




\' 


z-Koord. des Pinholes, d. h. axiale Auslenkung des PH 






Innendurchmesser des doppeltrichierfomiigen Pinholes 




5 


Offnungswinkel des Pinholes in Grad 




6 


Seitlicher, d. h. transversaler Neigungswinkel der Pinhole-Achse 


10 


7 


Axialer Neigurfgswinkel der PH-Achse 




8 


SchwachungskoeriizienL des Aperluniiatjerials in 1/ciu 





mpr . c 



Autoren: Dr. Nils Schramm ( FZ JUlich) 
Dr. Gemot Ebel (Fa. Scivis) 
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Datum: 20.08.2001 



MLEM-Variante der iterativen Multi-Pinhole-Rekonstruktion. Der 
Algorithmus beruecksichtigt die ortsabhaenglge Sensitivltaet und die 
ortsabhaengige Abbildungsf unktion einer MP-Apertur mit belieblg 
positionierten/geneigten Pinholes . 



20 



#include <stdiQ.h> 
# include <stdlib.h> 

include <sys/stat.h> 
ft include <unistd.h> 
# include <inatli.h> 
#include <-txnie.h> 

# in elude " / home / Schramm/ s r c/ ut i 1 / ut i 1 . h " 



30 



^define Wd2 M_PI 
#define Q 250 
#define U 2.5 

#define sci (m, k, 1) 
#define prj {1, k) 
#define rec{j,i,s) 
fidefine rot.(3,i,s) 
#define psf (h,l,k) 
#define kn2(l,k) 
fidefine kn3(l,k,h) 



/♦ Halber Winkelbereich */ 
/* Aufloesung der PSF */ 
/* PSF-Breite (98. "7%) */ 



sci t (m) -^-KL+tk) *-L+ll 
prj [ (i)*K4y:3 
rec[ ( j )'>-NS+(i)*S+s] 
rot I ( j )*NS+(i)*S+s] 
psf [ (h)*CC+(l+Cd2) *C+k+Cd2] 
kn2 [ ( 1+Ed2 ) *E+k+Ed2 ] 



/* Gemess. Projektionen 
/♦ Berechn. Projektion 
/* Rekonstruktion 
/•^- Rotierte Rekon. 
/* Abbildungsfkt. 

/* 2D-Kern 



kn3 [ ( 1+Ed2 ) *^EE+ ( k4-Ed2 ) *-E+h+Ed2] 



/* 3D-Kern 



Int N, S, M, K, L, NN, NS, KL, CC, Cd2; 

float P, A, p, W, cth, emu, gin, res; 

float *>:2, *y2, "^22, ♦hd, -^alf, ♦psi, *ther *amu; 



float • 
void 

void 

void 
void 

void 
void 
void 
void 
void 

void 

void 
void 
float 
void 



''PSF( 
SUM( 

PRJ( 

QUO{ 
REP( 

NRM( 
C0R( 
SAV( 
ROT( 
RAY( 

PIX ( 

SM2 ( 
SM3( 
LLH( 
DEV( 



float ^sen, 
float '*'prj, 
float *psf, 
float *prj, 
float *psf, 
float "^prj, 
float *rot. 



float 

float 
float 
float 
float 
float 



unsigned char *ise, 
unsigned char *ise, 

float *sen. 



*/ 
*/ 
*/ 
*/ 
*/ 



unsigned char *ise, unsigned char *ihb ); 
float *rot, float *sen, 
unsigned char *lhb ); 
float *rot, float *sen, 
unsigned char *ihb ) ; 
float ♦sci, int m ) ; 
float *ron, float ^prj , 
unsigned char *ise ) ;. 
float *cor, float *nrra ); 
*rec, float *cor ) ; 

*rec, float *rotf char *basen, int iter ); 
*rot, float *rec, float phi ); 

*V, float *H, float *a, float *b, float *c, float *c 
r, float t, float i2sico ) ; 

*F, float -^a, float *b, float ^-c, float *d, float r, 
float t, float i2sico, float ta, float ct ); 
float *sci, float *prj ); 
float *rec, float *rot ); 
float *prj, float ^sci, int. m ); 

FILE *devf, float *obj , float *rec, float llh, int n. ); 
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MAIN{): Hauptprogramm mit User-Eingaben, Speicherinitialisierung und 
der eigentlichen Rekonstruktionsschleif e . 
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mainO 
{ 

unsigned char *ise, '*^ihb; 

char scinl2563 , aptnt256], stanl256], objriI263, basenl2561, fulln[256]; 
int Sr ir j, n\, k, 1, t, stat, ticks, aflag, dflag, pflagr vflag, 1^ nsave; 
float *sci, *prj, ♦rec, *rot, -^cor, *nrin, *ron, *ser\, *psf, "^obj , llh; 
FILE *parf , -^aptf, *devf; 



==— 0) Bingaben vom Anwender 



printf( "XnProjektionen (float): " ); 
scanf( "%s", scin ); 

printf( ♦*Breite der Pro j . (pixj : " ); 
scanft "%i", &L ); 

printft "Schichten in Proj . [pix] : " ); 
scanf( "%i". ); 
printf( "Anzahl der Winkel : " ); 
scanf( "%i", &M ); 

printf( "Pixelgroesse in Proj. Imm] : " ); 
scanf( "%g"r &P ); 

pcintf( "Projektionen glaetten (nein«»0/ ja»l ) : •* ); 
20 scanf( "%i", &pflag ); 

printf( "Abstand RotAchse/Bildebene [mm]: " ); 
scanf( "%g"r iA ); 

printf( "Name d. Aperturmusters (text): " ); 
scanf( •'%s", aptn >? 
printf{ "\n" ); 

25 

printf( "Dicke des Kristalls Iinm] : " ); 
scanf( "%g", ficth ); 

printf( "TUDsorption d. Kristalls [l/cm]: " ); 
scanf ( "%g"r &cmu ); 

printf( "Intrinsische Aufloesung [iwa] : " ); 
30 scanf ( "%g"r &res ); 

print£( "Abs . Intensitaet d. Gammas "); 
scanf ( ••%g", &gin ); 
printf< -\n" ); 

printf( '^Startverteilung (float): " ); 
35 scanf ( *'%s", stan >; 

printf( "Breite der Rekon. Ipix] : " }; 
scanf ( "%i", ); 

printft "Schichten der ReV:on . Ipix): ); 
scanf ( "%i«, £S ); 

printf( "Voxelgroesse in Rekon. ivm] : " ); 
40 scanft ^^^g", &p ); 

printf( "Zahl der Iterationen: " ); 
scanf ( ); 

printf( "Volumen glaetten (n€in=0/ja=l ) : " ); 
scanft "%i'', ivflag ); 

printf( "Abbildungsf kt. voxelweise (nein=0/ ja-1) ?: 
45 scanf ( ••%i", fiaflag ); 

printf( "\n" ); 

printf( "Abvr. vom Original (nein-0/ ja=l) ?: " ); 
scanf ( tdflag ); 

if ( dflag 1 ) 
{ 

.printf ( "Originaldaten ( float ) : ) ; 
scanf ( "%s", objn >; 

) 

printf ( "Stammname d. Outputs: " ); 
scanf ( "%s", basen ); 
printft "Speichern alle?: " ); 
scanft "^i", finsave ); 
printft "Xn'' ); 



cis — rrr-n Ij Parameter wegschreiben ««»«=« 4^ i — = i ■ ■ ^ ■ ■ =as:««s ■■ . / 

sprintft fulln, "%s.par", basen ); 
parf « FOPENt fulln, "wt" ); 

fprintft parf/ "XnAlgorithmus: mpr.cSn" ); 

fprintft parf, "Parameterda-tei (text): %s.par\n", basen ); 
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fprintf< parf, "\n" 1; 

fprintfl parf, "Projektionen (float): %s\n", scin ); 

fprintf( parf, "Breite der Pro j . [pix) : %i\n", L ); 

fprintf( parf, "Schichten in Pro j . Ipix) : %i\n", K ); 

fprintf( parf, "Anzahl der Winkel: %i\n", M ); 

fprinttf{ parf, " Pixelgroesse in Proj . [nun]: %g\n", P ); 

^^fprintf( parf, "Pro j ektionen glaetten (nein=0/j a«=l ) : %i\n", pflag ); 

fprintf{ parf, "Abstand Bildebene/Rot^Achse Imra) : %g\n", A ); 

£printf{ p>arf, "Name d. Aperturmusters (text): %s\n*', aptn ); 

fprintft parf, "\n" ); 

fprintf( parf, "Dicke des Kristalls [mm): %g\n", cth ); 

fprintf( parf, "Absorption d. Kristalls [1/cm]: %g\n", emu ); 

fprintf( parf, "Intrinsische Aufloesung [mm}: %g\n'', res ); 

fprintf( parf, "Abs. Intensitaet d. Gammas [%] : %g\n", gin ); 

fprintf( parf, "\n" ); 



fprintf( parf, "Startverteilung (float): %s\n", stan ); 
fprintf( parf, "Breite der Rekon. Ipix): %i\n", N ); 

"Schichten der Rekon. [pix]: %i\n" , S ); 
"Voxelgroesse in Rekon. [mm]: %g\n'*, p ); 
fprintf( parf, "Zahl der Iterationen: %i\n", I ); 
fprintf( parf, "Volumen glaetten (nein»0/j a=l) : %i\n", vflag ); 
fprintf( parf, "Abbildungsfkt. voxelweise (nein«0/:ja*-l) ?: %i\n", aflag ); 
fprlntf( parf, "\n" ); 



fprintf ( parf, 
fprintft parf, 
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fprintf { parf, "Abw. vom Original (nein-0/ja«l)?: %i\n", dflag ); 

if ( dflag 1 ) 
{ 

fprintf ( parf, "Originaldaten (float): %s\n", objn ); 
fprintf ( parf, "Abweichungsdatei (text): %s,dev\n", basen ); 

} 

fprintf ( parf, "Rekonstruktion : %s.rnn\n", basen ); 
fprintf ( parf, "Spelchern alle?: %i\n", nsave ); 
fprintf ( parf, "\n" ); 



25 



30 



FCLOSE { parf ) ; 

/ = 2) Aperturmuster einlesen 

fprintf ( stderr, "Initializing. \n" )? 



aptf = FOPEN{ aptn, "rf ); 
fscanf( aptf, "%i", &T ); 

x2 MALLOCl T*sizeof ( float ) ) 
y2 « MALLOC( T-^sizeof ( float ) ) 
z2 « MAI,LOC( T*sizeof ( float > ) 
hd « MALLOC( T*siz€of ( float ) ) 
alf - MALLOC( T^sizeof( float ) 
psi = MALLOC( T*siseof( float ) 
the = MALt.OC( T*sizeof ( float ) 
amu - MALLOC( T*sizeof ( float ) 

for( t-0; t<T; t++ ) 
( 



Stat = fscanfC aptf, •%g%g%g%g%g%g%g%g", &x2ltl, ty2[t), fiz2[t], &hd[t3, 
&alf[t), &psi(t], &the[t], itamu[t3 ); 
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} 



if( Stat EOF ) 

ERROR( "MAINO: Fehlerhaftes Aperturmuster!' 



); 



alf It] « M_PI*alf [t)/180; 

psi It] -= M_PI*psitt)/180; 

theft] « M_PI'^theIt] /18C; 

amu(t] « amu[t]/10; 



/* Umrechnen in rad u. mm ^'Z 



55 



FCLOSE ( aptf ) ; 

/*= 3) Speicher anfordern und Felder initialisieren ass5s=sKSi==ss«=i=i=:'^ / 

if ( nsave < 1 I ( nsave > I ) ^ /*" nur sinnvolle Werte */ 
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nsave = I; 
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emu cmu/10; 
W ^ 2*Wd2; 
NN « N*N; 
NS N*S; 



/* Kristall'Mue in 1/mm */ 
/* Voller Winkelbereich */ 
/* ?Constanten *•/ 



sci ^ MALLOC( M*Fa.*si2eof ( float I ); 
prj « MALLOC< KL*sizeof{ float ) ); 
rec - MAIiI.OC( NN*S*siseof { float \ 
rot « MALLOC( NN-^S^sizeof ( float ) 
cor « MALLOC( NN^S^sizeof ( float ) 
nrm = MALLOC( NN*S*sizeof ( float ) 
ron « MALLOC( NW*S*siseof { float ) 
sen = MALLOC( T^Q^sizeof ( float ) ); 
Ise = MALLOC( T-^NN^S+sizeof ( unsigned char 1 ); 
ihb = MALLOC( T^NN^-S^sizeof ( unsigned char ) ); 

FREAD{ sci, sizeof( float >, M*i-CL, sctn ); 
FREADi rotr si2eof( float ), m*S, stan ); 

if ( pflag 1 ) 
SM2( sci, prj ); 

fort j<N; 3++ > 

for( i«0; i<N; i++ ) 
for( s=0; s<S; s++ ) 

recCj^NS+i^-S+s] « rot ls*NN+i*N+ j ] ; 

psf = PSF( sen, ise, ihb ); 

if ( dflag == 1 ) 

obj e MALLOCt NN*S*si2eof( float ) ); 

FREAD( rot, si£eof( float ), NN*S, objn ); 

for{ j=0; 3<N; j++ ) 
for( i=0; i<N; i++ ) 
for ( s-0; s<S: s++ ) 

obj [ j*NS+i^S + s] ^ rotls-^NN+i*N+j ] ; 

sprintft fulln, "%s.dev", basen ); 
devf - FOPEN( fulln, "wt" ); 



/* Speicher */ 



/* Einlesen */ 
/* Proj. glaetten */ 
/* reorganisieren ♦/ 

PSF berechnen 
/* Original lesen ♦/ 

/■*■ reorganisieren */ 
/* Abw.Datei */ 



»c»B=« 4 ) Rekonstruktion 
for( i=0; i<I; i++ 1 

printf( "Iteration %i von *i\n", i+1, I ); 



/* Schleife ue. Iter. 



ticks 



clockO; 
1 ) 



.if ( dflag = 
.. llh = 0; 



if( vflag == 1 ) 
SM3( rec, rot ) ; 

SET_FLT( cor, NN*-S, 0 ); 
SET_FLT( nrm, NN*S, 0 ); 

for( iwsO; ro<M; m++ ) 



/* Log-Likelihood */ 
/* Rekon, glaetten *^/ 
/* initlalisieren *•/ 
/♦ Schleife ue. Winkel */ 



{ 



fprintf< stderr, "%4i", m+1 ); 

SET_FLT( rot, NN*S, 0 ) ; 
ROT{ rot, rec, -m*W/M ); 

if { aflag 0 ) 

SUM( prj, rot, sen, ise, psf, ihb ); 
else 



/* Objekt ro tier en */ 
/* projizieren */ 
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PRJ{ pr^, ret, sen, ise, psf, ihb ); 











if < dflag «= 1 


> 

r»r^ iar*4 to. \ z 

^i.Jf J, f 111 J , 


/♦ Log'Likelh. */ 




• 








QUO( prj, sci. 


m ) ; 


/* Quotient ♦/ 


5 




m 






REP( rot, ron. 


prj , sen, ise ) ; 


/* rueckproj izieren */ 




* 






) 


ROT( cor, rot, 
ROT ( nrm, ron, 


m*W/M ); 


/♦ rueckrotieren */ 


10 








NRM ( cor , nrm ) ; 




Fa):t. normieren */ 










if 


[ vflag = 1 ) 
SM3( cor, rot 


); 


/■^ Fakt. glaetten */ 


15 








COR{ rec, cor ) ; 




/* Korrekturschritt */ 










print.f( "\nDauer: 


igXnXn",. (float) ( clock ( ) 


-ticks ) /CLOCKS_PER_SEC ) ; 










i£( dflag — 1 ) 

DSV( devf, obj, rec, llh, i+1 ); 


/* Abweichung */ 


20 








if 


( (i+l)%nsave = 
SAVi rec, rot. 


= 0 ) 
basen, 1+1 ); 


/* speichern */ 





SAV( rec, rot, basen, I ); 



letzten inuner ♦/ 



) 

PSF(): Berechnet die Sensitivitaet und die Punktbildbreite fuer das 
gesamte Rekonstruktionsvclumen. 

float *PSF( float *sen, unsigned char ''•ise, unsigned char "♦'ihb ) 
{ 

unsigned char *pis, *pih; 

int s, i, j, h, k, 1, t, Kmlr Lml; 

floai xO, yO, sO, xl, yl, y3, zB, X2, Z3, xlq, rdq, mrq; 
float Lmld2, KnxLd2, AdP, mdP, pnviP, oy, on, dy, dz; 
float de, fac, smi, sma, hmi , hma, hwb, *tse, ^pts, *pse; 
float yv, 2v, Yv, Zv, xvq, yvq, :-:r, yr, sr, br, xvr, yvr; 

float del, coa, cob, cog, bvq, bv, skp, nww, sum; 
float pina, pmq, n\pl, mlq, dq, hdq, req, *pef; 



Kmld2 = (?C-l)/2.0? 

Linld2 « (L-l)/2.0; 

Kml = K - 1; 

Ln\l « L - 1; 

xO = -p* (N-l)/2. 0; 
yO = -p^(N-l)/2.0; 
zO -p* {S-l)/2.0; 

req = res* res; 

mrq « p^p^-N^N/^ . 0; 

hmi = FLT^MAM; 
hma = -EXT_MA>:; 

for( t-0; t<T; t++ ) 
{ 

(A-xO)/(x2tt)-xO>; 
p*{inpl - 1); 
mpl*mpl; 



rapl 
pma 
miq 

pmq « pma* pma; 



dq = hd(t)*hd[t]^mlq; 
dq = dq'^cos( alf[tJ/2 ); 
dq =s dq + pmq + req; 
hwb sqrt( dq )/P; 



Konstanten ♦/ 



/* Untere Ecke des FOV */ 



/* Quadrat d. Intrins. */ 
/* Quadrat d. FOV-Radius ♦/ 

/* Min. und Max. der HWB */ 



/* Min/Max schaetzen */ 
Min. Vergroesserung */ 

/* Min. Punktbildbreite */ 



30 



35 



40 



45 



50 



60 



65 



3DOCID: <DE_ 



.10142421 Al_ 



DE 101 42 421 A 1 



15 



20 



25 



30 



35 



40 



45 



if ( hwb < hmi ) 
hrai ^ hwb ; 

mpl = <A+xO)/(:':2[tl+xO); 
pma - p* (mpl - D ; 
,mlq - mpl* mpl; 
pmq « pma* pma; 

xr = 1; 

yr = tan( psilt] ) ; 

sr = tan( t;he[tj ) ; 

br = sqrt( 1 + yr^-yr + sr*2r ); 

del acos( xr/br ); 

dq hdtt]*hdlt] *mlq; 
dq = dq + pmq; 

dq = dq/cos{ del + alfitl/2 ); 

dq - dq + req; 

hvrb = sqrt( dq )/P; 

if ( hwb > hma ) 
hma = hwb; 

> 

C = U-^hma; 

if ( C%2 0 > 

C++; 
CC « C*C; 
Cd2 = C/2; 

psf = MALLOCt Q*CC*si=eofl float ) ); 
for{ h'^O; h<Q; h++ ) 



/* Mas:. Vergroesserung */ 



Winl-iel PH/Kristall */ 



/* Max. Punktbildbreite */ 



/* Dimension d. PSF */ 



SpeicWer PSF */ 
/* PSF belegen */ 



sum = 0; 

hwb = hmi + h* (hma-hmi ) / tQ-1 > ; 

fori 1=-Cd2; l<=Cd2; 1++ ) 

for( k=-Cd2; k<=Cd2; k++ ) 

if ( 1*1 + k*k <«= U*U*hwb*hwb/4 ) 



/♦ max. 98.7% */ 



{ 



psf(h,l,h> exp< -2.772589* (1*1 + k*V:) / (hwb^-hwb) ); 
sum =« sum + psf(h,l,fc); 



> 

else 
I 

psf{h,l,k) ^ 0; 

) 

for( k=-Cd2; );<=Cd2; k++ ) 
for( l=-"Cd2; l<=Cd2; 1++ ) 

psf{hrl,k) = psf (h,l,k)/sum; 



/* normieren */ 



tse -= MALLOC< NN*S*sizeof I float > ); 



/* temp. Sen-Feld */ 



50 



55 



60 



pse = sen; 
pis «- ise; 
pih - ihb; 



/* schleife ue. Loecher */ 



for( t«0; t<T; t++ ) 

fprintf( stderr, "Calculating the PSF: %i/%i\r", t+1, T ); 

de = sqrt( hdlt}*(hdlt) + 2*tan( alflt}/2 )/amu[t]) ]; 
fac « 10000+gin*de*de/16; 
hdq « hdtt]*hdlt) ; 



oy « y2tt3/P + Lmld2; 
o2 - 22(tl/P + Knild2; 
dy = y2[t] - yO; 
d2 « z2It] - zO; 
AdP = (A - x2[t3 )/P; 
X2 •= x2lt3 ; 



/* vorbereiten */ 
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XV =s xo - X2; 
Yv « yO - y2Itl; 
Zv « zO - z2[t]; 

xr «= -1; 

yr e -tan( psi [t] ) ; 
zr = -tah^( the[t] ); 

= sqrt( 1 + yr*yr + sr*2r ); 
xr = xr/br; 
yr = yr/br; 
sr = zr/br; 

coa « cos( alftt)/2 ); 

>il « xO; 
pts « tse; 
smi = FLT_MA>:; 
sma = -FLT MAX; 



/*• Verbindungsvektor */ 



Richtung PH-Achse */ 



/*- Cos. d. PH-Oeffnung */ 
/*• Initiolisieren */ 



10 
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for{ j=0; j<N; j++ ) 
{ 

mpl = ( A-xl ) / (X2-X1 ) ; 
pma = p* { mpl - 1 ) ; 
mlq = mpl^'tnpi; 
pmq = pma'^pma; 
rodP = AdP/(X2-xl); 
y3 -" rodP^-dy + oy; 
23 « iiidP*'ds + o=; 
pindP = p*nidP; 

xlq - xl*-xl; 
yl -= yO; 



Proj, 



unteren Bcke 
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yv = Yv; 
xvq = xv*xv; 
xvr = XV* xr; 

for( i*0; i<N; 1++ ) 
{ 

1 = floor ( y3 ); 
23 « 23; 

rdq = xlq + yl*yl; 

2v — 2v? 
yvci » yv*yv; 
yvr = yv*yr; 

for( s«0; s<S; s+-i- ) 



{ 



if ( rdq < mrq ) 
{ 

- floor ( 23 ) ; 



innerhalb FOV? */ 
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if ( 1>«0 && KLml && k>«0 && k<Kml ) 
{ 

bvq «= xvq + yvq + 2v*2v; 
skp « xvr 4 yvr + zv*sr; 
bv = sqrt ( bvq ) ; 
cob «= skp/bv; 

if ( cob >« coa ) 



cog 
nww 



-XV /bv; 

1 - exp( -crau-^cth/cog ) ; 



/* in Prj? */ 



/* in PH? */ 

Wink. Lot */ 
/* Nachw.W. 
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*pts = fac*nvfw*cob/bvq; 

if ( *pts < smi > 

smi = *ptis; 
if( ♦pts > sma ) 

sma = "'■pts; 



/* Sensitiv. */ 
/* Min/Max d. Sen. ♦/ 
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) 

10 else 

{ 

*pts = 0; 
*pih « 255; 

) 

) 

15 else 
{ 

*pts «= O; 
*pih * 255; 

) 

) 

20 else 
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dq « hdq*mlq; /"^ I*och vergroessern */ 

dq = dq*cob; Winkel 2U Pinhole */ 

dq = dq + pmq; /* Voxelausdehnung V 

dq = dq/cog; Winkel zu Kristall V 

dq « dq + req; /* Intrins. Aufloesung */ 

hwb - sqrt:( dq )/P; /* Punktbildbreite */ 

♦pih = (Q-D* (hwb-hmi)/ (hma-hmi) + 0.5; 



{ 

♦pts « 0; 
*pih = 255; 

) 

::3 -= pmdP; 
zv 4e= p; 
pts++; 
pih++; 



yl 4- p; 
y3 pindP; 
yv +« p; 

> 

>a += p; 
XV p; 

) 

for( h^O; h<Q; h++ ) Sensitivitaet V 

*p3e - sini + h* (sma-sroi ) / (Q-1) ? 
pse++; 

) 

•fort i-0; i<NN*S; i++ ) /* Indizes d. Sen. 

45 < 

if( tseli) > 0 ) 

*pis = (Q-D* (tse[i)-sini) / (sma-smi) -v 0.5; 

else 

*pis = 255; 

50 pis++; 
} 

} 

printft -\n\n" ); 
55 FREEt tse )? 

return psf ; 

SUMO: Einfache Zentralstrahl-Projektion plus Faltung mit mittlerer 
60 Abbildungsfuntkion. Die Projektion beruecksi chtigt die ortsabhaer.gige 

Sensitivitaet. ^^^^^a,*.**.**.^***--^^*-^**-^**-*-*-/ 

void SUM( float *prj, float *rot, float *sen, unsigned char *-ise, 
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floaty *psf , unsigned char ♦ihb ) 



unsigned 'char *pls, *pih; 

int s; i, j, k, 1, t, kk, 11, IK, id>:; 

float xl, y3, z3, Z3; 

float Nmld2, Sn\ld2, Ln\ld2, Kmld2, AdP, pnvdP? 
float*DY, DZ, wl, w2, w3, w4 , tmp, mih, sum; 
float oy, o2, dy, dz, X2, ^tpr, *ppr, *pro; 
float *pse, *ppf, *Ppf; 

SET_FLT( prj, KL, 0 ); /* Pro j . nullsetzen */ 

tpr = MALLOC( ia-^si2eof( float ) ); /* temp. Projektion ♦/ 

Niftld2 - (N-l)/2.0; /* Konstanten */ 

Smld2 = (S-l)/2.0; 

?Cinld2 = (K~l)/2.0; 

Lmld2 = {L-l)/2:0; 

pis = ise; 
pih « ihb; 
pse sen; 

for( t=0; t<T; t++ ) /* Schleife ue . Loecher */ 

SET_FLT( tpr, Kl,, 0 )? Pro j . nullsetzen */ 

oy = y2 [ X ] / P -i- Lnild2 ; 

o2 = z2ltl/P + Kmld2; 

dy = y2[t)/p + Nmld2; 

dz = s2[t]/p + Snad2; 

AdP = (A - x2(tl )/P; 
X2 « x2It)/p; 

xl = -Nmld2; 
pro = rot; 

mih = 0; 
sum = 0; 

fort j-'O; j<N; j++ ) 



y3 = pmdP*dy + oy; 
Z3 = pmdP*ds + 03; 

for( i=0; i<N; i++ ) 
{ 

1 = floor { y3 ) ; 
DY = y3 - 1; 

z3 « 25; 

IK = 1*K; 

for( s=0; s<S; s++ ) 

tmp = ^pro; Objektwert * 

if ( *pis < 255 && tmp > 0 ) 
{ 



pmdP « AdP/(X2-xl)? 



/* Proj , d. unteren Ecke * 



tmp = tmp*pse [ *pis J ; 
k = floor ( z3 ) 7 
DS = tmp* ( z3 - k ) ; 



Sensitiv. dazu * 



w4 = DY*D2; 

w2 « tmp*DY - w4; 

vr3 - DS - w4 ; 

wl - tinp - DZ - w2; 



Gewichte 



ppr « tpr + IK + k; 
(•^ppr) +« wl; 



/*• auf sunwiieren * 
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ppr += 1; 
(*ppr) += w3; 

ppr +«= K; 

ppr — 1; 
(*ppr) += w2; 

mih ti«p*(*pih); E^SF mitteln */ 

10 sum +«* tmp; 

) 

23 -= pmdP; 
p»ro++ ; 
pi.s++; 

15 pih++; 

> 

y3 -= pindP; 

} 

20 xl++; 

idx = mih/sum +0.5; mittlerer Index */ 

ppf - psf + idx-CC; /* mittlere PSF */ 

ppr = tpr; 

for( 1=0; KL? 1+4 ) Faltung */ 

for( k-0; k<?C; k++ ) 
{ 

if I *ppr > 0 ) 
{ 

Ppf - ppf; 

fori Il«l-Cd2; Il<=l+Cd2; 11++ } 

fort hk=k-Cd2; kk<=k+Cd2; 1:1:++ ) 

^ if( •^Ppf>0 && 11>'=0 fc& IKL &i kk>«0 6fr kV:<K ) 
pr:J(llrkk) += (*ppr)*(*Epf ); 



25 



30 



35 



40 



45 



Ppf++; 



> 



) 

ppr++; 

} 

pse +« Q; 



free; tpr ); 

PRJO: Projektion mit voxelweiser Beruecksichtigung von Sensitivitaet 

void PRJ( float: *pr j , float *rotr float -^sen, unsigned char *ise, 
float *psf, unsigned char *ihb ) 

{ 

55 unsigned char **'pis, *pih; 

int s, i, j, \i, 1, l:k, 11; 

int 10. 11 r kO, kl, lOK, Lml , Vinl; 

float xl, y3, z3, 23; 

float Nmld2, Smld2, Lmld2, Kmld2r AdP, pmdP; 
float DYr DZ, wlr w2, w3, w4 , tmp; 
60 float oy, oz, dy, dz, X2, ♦pro, *ppf, *pse; 

float *ppO, *ppl, *pp2r *PP3, ^PP^; 

SET^FLT( prj, to-, 0 ); /* P^o^ • nullsetsen V 
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Nnad2 
Smld2 
Kmld2 
Linld2 
Kml ■ 
Lml 



(N-X)/2^.0 
(S-l)/2.0 
(K-l)/2.0 
(L-l)/2.0 
K - 1; 
L - I; 
* 

= ise; 
pih = ihb; 
pse «* sen; 



for( t=0; t<T; t++ ) 



/* Konstanten 



{ 



oy « y2Ct]/P + LmldZ; 
oz - z2lt)/P + tCmldZ; 
dy « y2lt]/p + Nmld2; 
dz « 22[t]/p + Snad2; 

AdP = (A - x2lt] )/P; 
X2 - x2[t3/p; 



/* Schleife ue. Loecher 



pro 



-Nmid2; 
rot; 



fori j«0; j<N; j++ ) 
{ 

pmdP « AdP/(X2-xl); 
y3 = pmdP*dy + oy; 
Z3 = pmdP^dz + os; 

for( i=0; i<N; i++ ) 
{ 

I = floor ( y3 J; 
DY = y3 - 1; 

10 « 1 - Cd2; 

II «= 1 + Cd2; 

s3 = Z3; 
lOK = 10*K; 



/* Proj . d. unteren Ecke 



for{ s-0; s<S; s++ ) 
{ 



•'pro; 



_if( *pis < 255 &£. tmp > 0 ) 



tmp = 


tmp*pse [*pis] 


n 


floor ( z3 ); 


D2 = 


trap-^CiS - k) ; 


kO = 


k - Cd2; 


kl = 


k + Cd2; 


w4 = 


DY*DZ; 


w2 = 




w3 « 


DZ - w4? 


wl = 


tmp - DZ ~ w2 ; 



Objektwert 



/* Sensitiv. dazu 



Gewichte 



/* PSF aufadd. 



ppf « psf + (*pih)*CC; 
ppO - prj + lOK + kO; 

for( 11^10; 11<«11; 11++ ) 
{ 

ppl « ppO; 
pp2 - ppO + K; 
pp3 = ppO + 1; 
pp4 * pp2 + 1; 



for( kk-kO; kk<-kl? kk++ ) 
{ 

if( *ppf > 0 ) 
{ 

if{ ll>-0 t& IKL && kk>*0 &£ kk<K ) 
♦ppl += wl* {*ppf ) ; 
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40 



45 



50 



if( Xl>«-1 K& IKLml 6& kJ:>«0 && kk<K ) 

*pp2 w2"^(*ppfl; 
if( 11>=0 && IKL && kk>=-l hk<haTa ) 

-^ppS +5= w3-^ (*ppf ) ; 
if( 11>=-1 && IKljnl && kk>«-l && ltk<Knil ) 

*pp4 +» w4*(*ppf); 

} 

ppf++; 
ppl++; 
pp2++; 
pp3++; 
pp4++; 

) 



ppO += PC; 

15 } 

} 

z3 -= pimdP; 
pro++ ; 
pis++; 

20 pih++; 

) 

y3 piitdP; 

} 

25 

) 

pse +« Q; 

} , 

) 

QUO(): Bildet den Quotienten aus gemessenen und berechneten Projek- 
tionen. pr j ( ) und sci { ) haben unterschiedliche Speicherorganisation. 



void QUO( float "'^prj, float *sci, int m ) 
{ 

int 1; 

for{ 1=0; KL; 1++ ) 
.for( k=0; k<K; k++ ) 
( 

if( sci(m,k,l] > 0 prjU/k) > 0 ) 

prjUrk) = sci(ro,krl)/pr3{l,k); 
else 

prj(l,k) « 1; 

) 

> 

REPO: Projiziert die fCorrekturvorschlaege in das (rotierte) Rekon- 
struktionsvolumen zurueck. Dies geschieht unter Beruecksichtigung der 
ortsabhaengigen Sensitivitaet. 

void REP( float *rot, float *ronr float *pr j ^ float *sen, 
unsigned char ♦ise ) 



unsigned char ^pis; 
int i, j, kr 1, t, IK; 
55 float xl, y3, sB, 23; 

float Nmld2r Smld2, Lmld2r Kmld2, AdP, pmdp; 

float DY, DZ, wl, w2, w3, w4, trap, quo; 

float oy, OS, dy, dz, 712, *pro, '^prn, *ppr, *pse; 

Nmld2 - (K--l)/2.0; Konstanten */ 

60 Smld2 = (S~l)/2.0; 

Kmld2 = (K-l)/2.0; 
Lmld2 = (L~a)/2.0; 



pis = ise; 
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for( i-0; i<NN*S; i++ ) 
if( rot[i) > 0 ) 

rot {13 - 0; 
else 

rottii] • -1; 

SET_rLTl ron, NN*S, 0 ); 

fort t=0; t:<T; t++ ) 
{ 

oy - y2rt]/P + Uiild2; 
oz - 22[t3/P + 1-0111032; 
dy - y2[t)/p + Nmid2; 
dz = z2[t]/p + SmldZ; 



AdP 
X2 



(A - x2tt) )/P; 
x2(tl/p; 



/♦ Aussenraun marhieren */ 



/* initialisieren -*/ 
/* Schleife ue. Loecher */ 



lU 



15 



xl = -Nnad2; 
pro •= rot; 
prn » ron; 

for{ j»0; j<N? j++ )' 
{ 

pmdP = AdP/(X2-xl); 
y3 « prodP-^dy + oy; 
Z3 « pmdP*ds + oz; 

for( i=0; i<N; i++ ) 
{ 

1 = floor ( y3 ) ; 
DY = y3 - 1; 



s3 
IK 



23; 
1*K; 



/* Proj . d. unteren Ecke */ 



20 



25 



30 



for( s=0; s<S; s+4 ) 
{ 

if( *pis < 255 &6 *pro != -1 ) 
{ 

k = floor t =3 ); 
^3 ~ k ; 



DZ 

w4 
w2 
w3 
wl 



DY^DZ; 
DY - w4; 
D2 - w4; 
1 - D2 - w2; 



) 



ppr - pr3 + IK + k; 
quo = wl* (*ppr) ; 

ppr +« 1; 

<guo w3* ( ^ppr ) ; 

ppr += K; 

quo +•=» w4*{*ppr); 

ppr 1; 

quo +« w2"*-(*ppr); 

tmp « pse t^'pis) ; 
* pro +— tw^* quo ; 
*prn trap; 



/♦ Gewichte */ 



/* mittlerer Quotient 



Sensitiv. */ 
Korr.Fakt. */ 
/"^ Normierung */ 



35 



40 



45 



50 



55 



z3 pmdP; 
pro++; 
prn++ ; 
pis++; 



60 
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y3 -« pmdP; 

} 

5 ) 

pse += Q; 

) 

) 

/ 4,+ ♦ *♦ 4.4- 4, + * 4.*. 4, * + * * + + ^ + + + ^. * ^4- ^ 

10 CORt): Ncrmierung der Korrekturf aktoren. 

void NRMt float *cor^ float ♦nrm ) 
{ 

int i; 

15 

for( i=0; i<NN*S; i++ ) 
if( nrmri] > 0 ) 

cor[i] ~ cor [i 3 /nrmli] ; 

} 

20 COR(): Anwendung der Korrekturf aktoren auf die aktuelle Rekonstruk- 

tion. 

void COR( float *rec, float ♦cor ) 
{ 

25 int i; 

for( i=0; i<NN*S; i++ ) 
if( recti] > 0 ) 

rec[i] = corIi3*-rec[i) ; 

) 

SAVn : Reorganisieren und speichern der Daten. 

void SAV{ float *rec, float *rot, char *basen, int iter ) 

35 char fulln(256]; 

int i, j, s; 

for( s=0; s<S; s++ > /* reorganisieren */ 

for( i=0? i<N; i++ ) 

for( j=0; j<N; j++ ) 

"~ rot (s*NN+i*N+j] « rec [ j*NS+i*S+s] ; 

sprintf( fulln, -*s . r%2 .2i" r basen. Iter ); /* Name -/ 

EWRITE( rot, si2eof( float ), NN*S, fulln ); /*• speichern */ 

ROTO: Dreht das Objekt ccw urn den Winkel phi. Die Objekte muessen 
inverse Speicherorganistion (j,i,sj besitzen. Merke: Das Zielvolumen 
wird in ROTO nicht initialisiert . 
»r**************++**4f ********+** *** + ********.***************^**.*^i.^* + ^t^ 

void ROT( float *rot, float *rec, float phi ) 

int i, j," 5, k, 1, ifac; 

float si, CO, i2sico, ta, ct, Nd2, Nd2iiil, trap; 
55 floats, y, r[9), at23r b[2], c[23, d[2], r, t; 

float AO, Al, BO, Bl, CO, CI, DO, Dl, fac, ox, oy, *pre; 
float *prl, *-pr2, ♦prS, *pr4, ♦prS, *pr6, *pr7, *-pr8, *pr9; 

fac = phi/(2*M_Pl); /♦ Auf 2^Fi reduzieren */ 

ifac = fac; 
60 fac = fac - ifac; 

phi = 2*fac*M_PI; 

if { phi < 0 ) /♦ Negative umrechnen */ 

phi -= 2*M_PI + phi; 
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if ( phi 0 ) * /* Vielfache von Pi/2 */ 

{ 

£or( j«0; j<N; j++ ) 
for( i-0; i<N; i++ ) 
for( s=0; s<S; s++ ) 
Lfi recCjrirS) > 0 > 
* rot<j,i,s) += rect3ri»s); 

* 

return; 

) 



ox = -co; 
oy " -si; 



if( M_F1 < phi && phi < 3*M_FI/2 ) 



10 



if ( phi == M_PI/2 > 
{ 

for( j«0; j<N; j++ ) 

for( i»0; i<N; i++ ) 

for< s«0; s<S; s++ ) 

if( rec(i,N-j-l,s) > 0 ) 15 
rot(j,i,s) rec(irK-j-lrS>; 

return; 

) 

if ( phi =« M_PI ) 20 
{ 

• for( j-0; j<N; j++ ) 

for{ i«0; i<N; i++ ) 
for( s=0; s<S; s4-+ ) 

if( rec{N-j-l,N-i-l,s) > 0 ) 

rot(j,i,s) +« rec<N-j-l,N-i-lrS); 

return; 

) 

if( phi = 3*M_PI/2 ) 
{ 

for( j=0; j<N; j-f-f ) 

for< i=0; i<K; i++ ) ... 
for( s«0; s<S; s++ ) 

if( rec(N-i-l,j, s) > 0 ) 

rot(j,i,s) += rec(N-i-l, j,s); 

ret:urn; 

) 



25 



30 



35 



Nd2 = N/2.a; 
Nd2ml - Nd2 - 1; 
si = sin{ phi ); 
CO = cos { phi ) ; 

if ( 0 < phi && phi < M_PI/2 ) /* Welcher Quadrant? */ 

{ 

AO « -Nd2+<co - si); 
Al « -Nd2* (si + CO); 

ox = -si; 
oy « co; 

) - - 

if ( M_FI/2 < phi && phi < M_PI ) 
{ 

AO = -Nd2*co + Nd2ml*si; 
Al » -Nd2*si - Nd2inl*co; 

phi « phi - M_PI/2; 55 
si * sin{ phi ) ; 
CO « cos { p»hi ) ; 
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AO « -Nd2na*(co - si ) ; 
Al « -Nd2ml*(si + co); 

,phi - phi - M_PI; 
5 si « sin ( phi ) ; 

CO = cos{ phi ); 

ox = si; 
oy « -co; 

} 

10 

if ( 3*M_PI/2 < phi Afi phi < 2*M_PI ) 
{ 

AO = -Nd2ml*co + Nd2+si; 
Al - -Nd2ml*si - Nd2*co; 

15 phi = phi - 3*M_PI/2? 

si « sin( phi ); 
CO - cos < phi ) ; 

ox = co; 
oy = si; 

20 ) 

BO AO - si; 
Bl -= Al + co; 
CO « BO + co; 
CI ■« Bl + si; 
25 DO = CO + si; 
Dl « CI - co; 

i2sico = l/(2*si*co); 
ta = si/co; 
ct « 1/ta; 

30 

pre *= rec; 

for( j»0; j<N; j++ ) /* Objekt rotieren */ 

{ 



35 



40 



45 



a[0] 


EZ 


AO; 


a[l] 




Al; 


b[01 




BO; 


bfll 




Bl; 


ctOJ 




CO; 


c[l] 




CI; 


dIO] 




DO; 


dtl] 


St 


Dl; 


for ( 


i«0; i<N; i4-+ ) 


{ 






X 




0.5*( a[0) + c{03 }; 


y 




0.5*( a[l) + c[lj ); 


r 




floor ( y ) ; 


t 




floor t X ) ; 


k 




r + Nd2; 


1 




t + Nd2; 



50 pr5 « rot. + 1*NS + k*'S} 

PIX( F, a, h, c, d, r, t, i2sico, ta, ct >; 

for( s«0; s<S; s++ ) 
{ 

55 trap " *pre; 

if ( tmp > 0 ) 

{ 

pr2 = pr5 - S; 
pr8 - pr5 + S; 
60 pr4 = pr5 - NS; 

prl = pzA - S; 
pr7 = pr4 + S; 
pr6 = pr5 + MS; 
pr3 = pr6 - S; 
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p;:9 « pr^ + S; 

*prl += FtOl^tn^? 
*pr2 +• r{l)*tinp; 
♦pr3 ri23^tinp; 
*pr« r[3)*tn^; 
* -^prS += ri4]*tmp; 

■'-pr? +» ?|6]*trap; 
*pr8 4« r[7)*tmp; 
*pr9 4- Fiej*tmp; 

) 

pre++; 
pr5++; 



atO) 


+= 


ox; 


a[l] 


+= 


oy; 


b{0] 




ox; 


b[l] 


+= 


oy; 


cro] 


+= 


ox; 


c[l] 




oy; 


dl03 




ox; 


d[ll 




oy; 



AO += oy; 
Al -= ox; 
BO +- oy; 
Bl -« ox; 
CO +« oy; 
CI -= ox; 
DO +■» oy; 
Dl -« ox; 

> 

} 

RAY{); Berechnet die Schnittf laechen swischen einem rotierten Pixel 
und den drel ueberstreichenden Proiektionsstrahlen. 

void RAY( float *V, float float ^a, float *b, float *c, float ♦d, 

float r, float t, float i2sico ) 

{ 

float dl, d2r d3r rpl, tpl; 

rpl = r + 1; 
tpl « t + 1; 

ifi c[ll > rpl ) 
{ 

d2 - c(13 - rpl; 

if ( bin > rpl ) 
{ 

dl - btl) - rpl; 

V12) « <d2*'d2-dl*dl)*i2sico; 

} 

else 
< 

if( d{l} > rpl ) 
{ 

d3 = d[l] - rpl; 

Vt2] « Id2*'d2-d3*d3) *-i2sico; 

\ 

else 
t 

V[2] d2*d2*i2sico; 

} 

) 

) 

else 
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Vt2] - 0; 

} 

if{ atH < r ) 
5 i 

d2 « r - all) ; 

if( dll] < r ) 
{ 

dl - r - dtl); 
10 V[03 = (d2*d2-dl*dl)*12sicc; 

) 

else 
{ 

if( t>[l] < r > 
{ 

15 d3 = r - bill; 

V[0] ^ (d2*d2-d3*d3) -^123100; 

) 

else 

V[03 * d2*d2*i2sico; 

20 ) 

) 

) , 

else 

{ 

VIO) = 0; 

25 ) 

Vtl] « 1 - V[01 - V[2); 

if ( blO] < t ) 

30 d2 t - b[01; 

.if{ aid < t ) 
{ 

dl = t - alO]; 

H[0] = (d2*d2-dl-^dl)*i2sico; 

} 

else 

{ 

if I c[01 < t ) 
{ 

d3 t - ctO); 

HlOl « (d2*d2-d3*d3)*t2sico; 
} 

else 
{ 

H(0] « d2^d2*i2sico; 

) 

else 
{ 

H[03 0; 

} 

if ( d[01 > tpl 1 
{ • 

d2 « d[0] - tpl; 

55 if( clO) > tpl ) 

< 

dl ctO] - tpl; 

H[21 « {d2*d2-dl*dl)*i2sico? 

} 

else 
60 { 

if( a[0) > tpl ) 

{ 

d3 = a[0] - tpl; 

Ht2] « (d2'^d2-d3*d3)*i2sico; 
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) 

else 
{ 

HI2J = d2*d2*i2sico; 

) 5 

) • 

else 
{ 

Ht2] « 0; 

} l« 

HID = 1 - H[03 - H12J; 

) 

PIX(): Berechnet die Schnittf laechen zwischen dem rotierten Pixel 

und den neun unterliegenden Pi>:eln eines nicht-rotierten Rasters. 15 

void PIX( float ^F, float *a, float ♦b, float *c, float *d, float r, 
float t, float i2sico, float ta, float ct ) 

float dx, dyr v, h, v[3), Ht3], tplr rpl; 

tpl « t + 1; 
rpl «= r + 1; 

F{0] - F[2] - F[6] = F[81 - 0; ^. 
RAy( V, H, a, b, c, d, r, t, i2sico ); 

if( a[0) t &S. all3 <= r ) / Fall 1 und 2 */ 

{ 

dx «= t - alOl; 
dy » r - atl); 
V ss dy - dx^ta; 
if( V > 0 ) 
{ 

h cJx 4 dyta; 
F[0) « 0.5*( v*dx + h*dy ); 

1 

else 

FlO] - dy-dy*i2sico; 

} 

if ( b(0) <« t && bll] <= r ) Fall 3 und 4 */ 4q 

{ 

d>: - t - btO); 
dy - r - b[13; 
h * dx - dy*ct; 
if ( h > 0 ) 

( 45 

V « dy + dx*ct; 

F[0] = 0.5* ( v*dx + h*dy ); 

} 

else 

FfO] = dx*dx*i2sico; 

) 50 

if( a(0] > t &i atl) r && btO] <= t && b[13 > r ) /♦ Fall 5 */ 

{ 

dx « a[0} - t; 
dy « r - a [1) ; 

V •= dy - dx*ct; 55 
if { V > 0 ) 

F[0) - 0.5*v*v*ta; 

} 



if I b[0) t &i b(l) > rpl ) Obere, linKe Ec);e */ 

( 

.dx « t - blOl ; 
dy « bCD - rpl; 
h » d>: - dy*ta; 
if ( h > 0 ) 
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{ • 

V = dy + dx*ta? 

F[6} « 0.5* t v*dx + h*dy 1; 

> 

S else 

F[61 ^ dx*dx*i2sico; 

) 

ifl ctO) <« t && c[l] > rpl ) 
{ 

It) = t - cJO] ; 

dy «= c[l} - rpl; 
V B dy - dx*ct; 
if ( V > 0 ) 
{ 

h = dx + dy^ct; 
15 F[6] = 0.5*( v*dx + h-^dy )? 

} 

else 

F[6l = dy'^dy*i2sico; 

} 

20 iff bl03 t && b{13 <= rpl && c[0) > t &i c[l] > rpl ) 

dx » t - btO); 
dy « rpl ' b(13 ; 
h - dx - dy*ct; 
if( h > 0 ) 

F(61 « 0.5*h*h*ta; 

) 

if( clOl > tpl &c cri) > rpl } Obere, rechte Ecke */ 

{ 

dx « c[Oj - tpl; 
dy - cllj - rpl; 
V = dy - d::*ta; 
if I V > 0 ) 
i 

h d>: + dy*ta; 
F(e] « 0.5- ( v*dx + h*dy ); 

) 

else 

F{S1 = dy*dy*i2sico; 

} 

if { dtO] > tpl &i d(ll > rpl ) 
{ 

dx - dlOl - tpl; 
dy « dtl] - rpl; 
h = dx - dy*ct; 
if ( h > 0 ) 
{ 

V = dy + dx-^ct; 
F[3) = 0.5* ( v+dx + h+cy ); 

} 

else 

F{8J » dx*dx*i2sico; 

50 ) • 

if { clO] tpl && clD > rpl && dtOl > tpl && d[l) <»= rpl ) 

{ 

dx « tpl - cEO); 
dy « c[l] - rpl? 
55 V = dy - dM"^ct; 

if ( V > 0 ) 

F(8) « 0.5*v*'v*ta? 

} 

if ( d[0) > tpl && d[l] r ) /* Untere, rechte Ecke */ 

60 ( 

dx * dtOl - tpl; 
dy « r - d[ll; 
h « dx - dy+ta; 
if ( b > 0 ) 
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{ 

V -* dy + d2e*ta; 

F12) 0.5*( v*dx + h*dy ); 

) 

else 

FI2) « dx*dj«*-i2sico; 

} * 



if( alO] > tpl && all] r ) 
{ 

dx a [0} - tpl; 
dy = r - atl] ; 
V dy dx*ct; 
lf( V > O ) 
{ 

h « ds + dy*ct; 

Ft2} « 0.5* { v*dx + h*dy ); 

) 

else 

F12] = dy*dy*i2sico; 

} 



10 



15 



if( al03 <« tpl atl] <= r &fi d[0] > tpl 6& dll] > r ) 
{ 

dx = d|01 - tpl; 
dy = dill - r; 
h » dx - dy*ct; 
if< h > 0 ) 

F12] « 0.5*h*h*ta; 

> 



20 



Ftl] = V[03 - r[0] - F[2) 
F(7] = V[2] - F16] - F[8) 
F(3) - HIO) - F[0] - F16] 
FC51 = H123 - F[23 - FI83 
F[43 » 1 - HIO) - HI2] - 



30 



SM2(): 2D-Glaettung (binomial) der gemessenen Pxoj ektionen. 



void SM2 ( float *sci, float *pr:j ) 
{ 

int m, h, 1, kk, 11, Z, Ed2; 
float *kn2, sum; 

E =^ 3; 
Ed2 = E/2; 



35 
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kn2 = MALLOC( E*E*si2eof( float ) ); 



/* Kern */ 



kn2(-l,-l) = kn2(-l,l) «= kn2{l,-l) = kn2(l,l) = 1; 
kn2(-l,0) - kn2<0,-l) = kn2l0,l) = kn2(l,0) - 2; 
kn2(0,0) « 4; 



45 



for( xn«0; «)<M; m++ ) 
{ ' 

for( 1=0; KL; 1++ ) 
fort k=0; k<K; k++ ) 
{ 

sum = 0; 

prj (l,k) « O; 



/* Glaettung */ 



50 



for( hk=-Ed2; kk<=Ed2; kk++ ) 

for( 11 — Ed2; ll<»Ed2; 11++ ) 

if( 1+11>«0 && l+ll<i. fifi 3:+kl:>=:'0 &£ k+kk<K ) 
{ 

sum kn2{ll,kk); 

prj II, k) += kn2(ll,kk)*sci(itirk+kk,l+ll); 

} 



55 
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prj{l,k) = prj (1, k) /sum; 
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fori 1=0. KL; 1.. ) ^* uebertragen V 

for{ \i'=0; k<K; k++ ) 

sci(m,k,l) « prjdrk); 

) 

FREEt kn2 ); 

SM3(): 3D^Glaettun.^(binomiall^ 



II) 



/♦^ Kern */ 



/* Glaettung */ 



void SM3( float *rec, float *-rot ) 

^ int s, i, jr k, 1, E, Ecl2, EE; 
float *kn3, sum; 

15 E » 3? 

Ed2 « E/2; 
EE == E*E; 

kn3 - MALLOC( EE-^E^siseof ( float ) ); 
SET_FLT( kn3, EE*E, 0 ); 

20 

kn3(0,0r0) « 2; 
y:n3(0,0,-l) = 1; 
kn3t0,0,l) « 1; 
kn3l0,-l,0) - 1; 
kn3(0,l,0) o 1; 
-5 kn3(-l,0,0) = 1; 

kn3(lr0,0) = 1; 

toti j=0; j<N; j++ > 
for( i=0; i<N; i+4 ) 
for( s=0; s<S; s++ ) 
30 if( rec(j,irS) > 0 ) 

{ 

sum = 0; 
rot(j,irSl = O; 

tori 1«-Ed2; l<«Ed2; 1++ ) 
35 k«-Bd2? k<=Ed2; k++ ) 

i4-k<N && s+h>«0 && s+h<S ) 

sum += kn3(lfkfh>; 

40 rot(j,irS) 4« kn3{lrk.h)*rec(j+l,i+k,s+h); 

} 

rot<j,irS) - rott3,i/S)/suni; 

> 

45 , /♦ uebertragen */ 

fori i=0; i<NN*S; i++ ) , ^ 

recti] = rotlij; 

FREE! Vcn3 ); 



} 

50 / 



/ 

float LLH( float *pr j , float ♦sci, int m ) 
{ 

55 int kr 1; 

float llh, tmpl, tmp2; 

iin = 0; 

for( 1=0; !<!>; 1++ ) 
60 for( k-0; k<K; k++ ) 

* tmpl » prj (Irk); 
tmp2 - sci(m, krl); 
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if ( tnipl > 0 i& i:inp2 > 0 ) 

llh +« tmpa-^log (tmpl ) - tmpl; 

} 

return llh; 5 
K... 

* DEV(): Berechnet die mittlere quadratische Abweichung zwischen ObjeV;t 
und Rekonstruktion. ^ 

LO 

void DEV( FILE ♦devf. float *obj, float *recr float llh, int i ) 
{ 

int \i, cnt; 
float dev, tmp; 

cnt « 0; 15 
dev — O; 

for( k«0; k<NN*S; k++ ) 
if ( rec[k) > 0 ) 
{ 

cnt++ ; 20 
tmp » obj tk] - rec[k); 
dev += tmp* tmp; 

} 

dev = dev/ cnt; 

25 

fprintf{ devf, "%71%15.6g%15.6g\n", i, dev, llh )? 
fflush( devf ); 

30 



util.h Autor: Nils Schramm 35 

Datum: 10.12.1998 

Headerdatei zu util.c. 

#ifndef _-«til 40 

ifdefine util 

^include <stdio.h> 
^include <stddef.h> 
include <stdlib.h> 

#include <sys/stat.h> 45 
g include <unlstd.h> 
# include <raath.h> 
# include <limits.h> 
g include <float.h> 

fidefine TRUE 1 50 
♦define FALSE O 

void ERROR ( char *text ); 

FILE *FOPEN( char *name, char *nK>de ); 

void REWIND { FILE -^file ); 

void FCLOSE( FILE ♦file ); 

int . FLEN_BIN« FILE *file ); 

void FREADl void -^buffer, int size, int num, char *name ); 

void FWRITE{ void "*^buffer, int size, int num, char *name ); 

void FAPPEND( void *buffer, int size, int num, char *name ); 

void *MALLOC( int bytes ); 

void FREE ( void *ptr ) ; 

int MAX__INT( int *data, int num ); 

int MIN_INT( int *data, int num ); 

float MAX_FLT( float *data, int num )^ ^ 
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float MIN_FLT( float Mata, int num ); 
int SUM_INTt int *data, int num }; 
float SUM FLT{ float *data, int num ) ; 
void SET'iNTC int *data. int num, ^^^^^l >; 
void , SET"fLT( float *data, int num, float val ) , 
void ' COPY INT< int *dest, int -sour, int num ), 
«^srt cOPY~FLT( float *dest, float ^sour, xnt num 

void SUB_FLT( float "'dest, float souxj., iaoou 

int FLEN_ASC( FILE "-file ); ^ 

void READ INT ( int *data, int nslx, nrow, ^^^^^^^^ .,,3^^ j. 

void WRITE,INT( int -data, int nsli, nrow ^f^^f^^^^^^^ ); 

void APPEND_INT{ int *data, xnt nsli, xnt nrov xnt n ^^^^^ 

void READ_rLT( float -data, ^^t nsli, ^^^t^^ow, int ncoL, char *name ); 

void WRITE FLT{ float *data, int nsii, ^^^J^^^^' ^^nt ncoL char *nanie ); 

void APPEND_FLT( float *dat:a, int nsli, xnt nrow, int ncoi. cn 

int IS_BIG_ENDIAN ( void ); 

void SWAP TWO( char *chrl, char *chr-- ); 

void SWAPT void *ptr, int size ); 

20 ffendif 

— *v 



10 



15 



— * — rt:;r;rirrc;;:r* 

25 util.c Datum: 10. 12. 1998 
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^include "util.h" 



ERRORO : Standard Fehlerbehandlung^ ^ ^^^^^^^^^ 

void ERROR ( char "^text ) 

fprintft stderr, "\n%s\n", text ); , . „ . 
fprintf{ stderr, "E>:iting to system. .. \n\n 
exit( 0 ); 

FOPENt): Oeffnen einer ^^^^f J;;^^^^^^^,,^^^^.****^**^^**-*.****^**^****/ 

FILE *FOPEN( char *naroe, char *mode ) 
{ 

FILE •'file; 

file = fopen( name, mode ); 
if ( file NULL ) 

ERROR ( "FOFENO: Can't open file*' )? 

return file; 

REWIND( ) : Zurueckspulen ^^^^^^^f ^^i;^^,^^^^,^^^**^**^.***^** *******^'*/ 

void REWIND { FILE *file ) 
{ 

rewind ( file ) ; 

— * — 

60 rc--OSE(U Schltessen einer 

void FCLOSB( FILE *file ) 
C 
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f close ( f ile^ ) ; 

) 

rLEN_6IN(): Bestlmme die Laenge einer Date! in Byties. 

int FLEN»_BIN( FILE *file ) 
{* 

int len; 

struct Stat fst; 

£stat( £lleno( file ), &fst ); 
len = fst. St size; 



29 
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return len; 

} 

FRE:aD(): Blnaeres Lesen von Daten. 

void FREAD( void *bufferr int size, int num, char *name ) 
{ 

int len; 20 
FILE *fil€; 

file = FOPENl name, "rb" ); 
leri = FI.EN_BIN( file ); 

if( nuin*si2e len ) 25 
ERROR I TREADO: Wrong file size" ); 

£read{ buffer, size, nuin, file ); 

FCLOSE( file ); 

) 30 
FWRITE(): Binaeres Schreiben von Daten. Datei wird neu erzeugt, 

void FWRITE{ void *buffer, int size, int num, char *name ) 
{ 

FILE *file; 

file « FOPENl name, "wb" ); 
fwrite { buffer, size, num, file ) ; 
FCLOSE( file ); 

) 

y 4. 4p 4r «. ^ ^ 4. 4. 4^ 4> ••- «■ «■ "fk 4' «> 4- « <* 4. 4. ♦ «• 4^ -f- « « 

FAPPENDO : Binaeres Schreiben von Daten. Daten werden an eine be- 
stehende Datei angehaengt. 

void EAPPEND( void *buffer, int size, int num, char *naine ) 
FILE * file- 
file - FOPEN{ name, "ab" ); 
fwrite ( buffer, size, num, file ); 
FCLOSE( file ); 

> 55 
MALLOC(): Anfordern von Speicher. 

void *MALLOC( int bytes ) 

{ 60 
void *ptr; 

ptr = malloc( bytes ); 
if ( ptr NULL 1 
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ERROR! "MALLOC-O: Can^t allocate memory" ); 
return ptr; 

, ) * * 

-^^^ — — / 

void FREE( void *ptr ) 
( 

10 free( ptr ); 

^ MAX into: BestiKmte das Maximum eines Integer^ 
********'^** 

15 int M7^^ INT( int *data, int num \ 
{ 

int i; 
int mayii 



20 



max - INT^MIN; 

for( 1=0; Knurn; i++ ) 
if ( data[i] > max ) 
max data [11; 

return mas;; 
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int MIN^INTl int *data, int num ) 
{ 

int i; 
int iriin; 

min « IMT_MAX; 

for( 1=0; Knum; i++ ) 
if ( data(ll < min ) 
min = datati) i 

return min; 



MAX.rLTO: Bestinane das Maximu™ eines Float-^^ ^ 

float MAX_FLTl float *data, int num ) 
{ 

int 1; 
float max; 

max « -FLT^MAX; 

fort i**0; Knum; i++ ) 
ift data[ij > ma?: ) 
max « datati] ; 

return max; 

*****•'***•****■*•***********''***** 

float MIN_FLTt float ^data, int num ) 

/ 

60 int 1; 

float min; 

inin - FLTJW:; 
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for( i«»0; Kfium; i++ ) 
if( dataCiJ < min ) 
min « data [LI ; 

return min; 

^ SUM_INT(): Summiert einen Integer-Datensatz, 

int SUM_INT( int *data, int num ) 10 
{ 

int i; 
int sum; 



sum 0; 

for( i=0; i<num; i++ ) 
sum «= sum + data[i]; 
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return sum; 

} 

+ + * + + * *'V** + **.*4.<*r*** if 20 

SUM_FLT(): Summiert einen Float-Datensats. 
float SUM_FLT( float *data, int num ) 

^ int 1; 25 
float sum; 

sum « 0; 

for( i=0; Knum; i++ ) 
sum » sum datafi]; 

return sum; 

) 

y 4. + + + 4- + + -J^ + * + 4- + + .^r -ir-*- * 

SET_INT{): Initialisieren eines Integer-Datensats. 

void SET__IMT< int *data, int num, int val ) 
{ 

int i; 

for( i=0; Knum; i++ ) 
.d3ta[i] - val; 

) 

/<<-*-<^+ + + 4r4tV**+**>4.*<^4-*4r*4.4'4.4-4'4'+*- + *4-4-**4r*4'*4-+*-4-+4.*4.^4'+ ■♦- + 4'4' + 4:'*'4'4- 

SET_FLT(}: Initialisieren eines Float-Datensatz. 

void SET_FLT{ float *data, int num, float val ) 
{ 

int i; 

for( i«0; Knum; i++ ) 50 
datali] val; 

) 

/4^4-^4r4k-4r4.<^^4rA><fc^«4'<*'«4'4>4'>^<»<^4.4p^*4f4F4>4k>4rfr<4>4'<V4r<lt>4^4'4'4r4>4'4-4.*4>4>4.«. 

COPY_INT{): ?Copieren eines Integer-Datensatz. 

4'4-4'4-'*-4.4.<<.4-4'4-4^4-4-4'4^<t-4-4-<A-4--fr4.^^4.4.4--^4-4.4-4.4r4>4'4..^4-f-4.4--V4-4.4'4-«^« / 
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void COPY_INT( int *dest, int *sour, int num ) 
{ 

int i; 

for( i»0; i<num; i++ ) 

dest[il =» sour[i]; 60 

) 

COPY_FLT(): Kopieren einens Float-Datensats. 

♦ 4-4.*.4.4.4.4-**»*4.**4f-^+-**4.«.4.4r4r*''*-****4.4.*4.4.4.4.**4.*4r'*4.* / 
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void COPY_FLT( float *dest, float *sour, int nuin } 
I 

tnt i; 

fori i«0; Knum; i++ ) 
destU) - sourli); 



10 



15 



ADD INTO: Addiert swei I^teger-Datensaetze^ 
void ADD^INTl int *dest, int -sourl, int *sour2, int nxm ) 



int i; 



fori i«0; i<num; i++ ) 

destli] - sourltil + sour2li]; 



...^^iiiJ:^^^'^^^ ^ * 

20 void ADD^n.T( .loat *dest, float ^souri. float *30ur2, int nu. ) 



{ 

int i; 

fort i«=0; i<num; i++ ) 

destU3 = sourltil + sour2[il; 

SUB INTO: Subtrahiert ^^.ei Integer-Datensae^^^ 

void SUB^INTt int -dest, int ^sourl, int *sour2, int num ) 
30 { 

int i; 



for{ i»=0; i<nuni; i++ ) 

destlll « sourltil - sour2li3J 



35 



40 



' * * ****** 



void SUB_FLT( float *dest, float *sourl, float ^sour2, int num ) 



{ 

int i; 



fori x=0; i<num; i++ ) 

destlil « sourltil - sour2[il; 



50 
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int FLEN_ASC( FILE *file ) 
{ 

int status, cnt; 
float tmp; 

cnt « 0; 

while ( 1 ) 

^ status = fscanf{ file, "%g"r ^tmp >; 
if( status EOF 1 
break; 



cnt++f 

60 > " 



REWIND ( file ); 
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return cnt; ^ 

) 

READ_INT{): Elnlesen eines ASCXI-Datensatizes in elne Integer-Matrix. 
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ . 

void,^RErfb_INT ( int ^data, Int nsli, Int nrow, int ncol, char ^name ) 
{* 

int 1, len; 
float tmp; 

FILE *file; 1q 

file ^ FOPENt name, "rt" ); 

len «= FLEN_ASC( file ); 

if ( nsli-^nrow-'-ncol != len ) 

ERROR( "READ_INT(): Wrong file size" ); 15 

, for( i=0; i<nsli-*nrow*ncol; i++ > 
{ 

f scant ( file, "%g*'r fitmp ); 
,if ( Ziap >« 0 ) 

data[i) = (int)tmp + 0.5; 20 
else 

data[i) « (int) tmp - 0.5; 

) 

FCLOSE( file ); 

} 25 

WRITE_INT(): Schrelben eines Integer-Datensatzes im ASCII-Format. 
Date! wird neu erzeugt. 



void WRITE_INTl int ^data, int nsli, int nrow, int ncol, char *naroe ) 
{ 

int i, j, k; 
FILE *file; 

file « FOPEN( name, "wt" ); 

for( k-0; k<nsli; k++ ) 
K 

for( i-0; i<nrow; i++ ) 
{ 

for( j=*0; Kncol; ) 
( 

fprintfl file, "%i\t", data [k*nrow^ncol+i*ncol+j 1 ); 

. 1 

fprintf( file, "\n" ); 

) 

fprintfC file, "Xn" ); 

) 



FCLOSE( file ); 
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APPEND_INT ( ) : Schreiben eines Integer-Datensatzes im ASCII- Format. SO 
Daten werden an bestehende Datei angehaengt. 

void APPEMD_INT( int *data, int nsli, int nrow, int ncol, char *name ) 
{ 

int i, k; 55 
FILE -^file- 
file = FOPEN( name, "at" ); 
for( k=0; k<nsli; k++ ) 

{ .60 
for( 1«0; i<nrow; i++ ) 
{ 

for( j-0; j<ncol; j++ ) 
{ 
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fprintff file, "%i\t". data[k*nrow*ncol+i*ncol+1J ): 
£printf( file, "\n ). 
fprintfl file, "\n" ) ; 



FCLOSE( file ); 
,0 ^**;;;;**^;or.lnle3eneine3^.SCll;^^.en^ 

void REAO.rx,T( float *data, int nsli. int nrow, int ncol, char *na»e , 

{ 

int i, len; 
15 float tmp; 

FILE * filer- 
file - FOPEN( name, "rt" ); 

ler. « FLEN_KSC{ file ); 
20 if( nsli*nrow*ncol !- len ) ^.^^ 

ERROR{ "READ^FLTO: Wrong fxle size 

for( i*Of i<nsli*nrow*ncol; i++ ) 

^ fscanf( file, &tinp ); 

25 data I i] = 

) 

FCLOSEC file ); 

) 
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void «RITE_n.T( float *data, int nsli, int nrow, int ncol, char *na»e , 

* int ir jr H; 
FILE *file; 

file « FOPEN( name, "wt" )? 
for( k=0; >:<nsli; k++ ) 

for{ L«0? Knrow; i++ ) 

^ fort j=0; j<ncol? 1 ) 

^ fprintf( file, "%g\t", dataU-nrow-ncol+i^ncol+j) ); 

fprintfl file, "\n" ) ; 

fprintf* file, "\n" ); 

J - - 

FCLOSE( file ); 

oaten werden an bestehende Datei angehaengt.^^^^^^^^^^^^^^^^^^^^,,„/ 
void APPEND.n-T, float *data, int nsli, int nrow, int ncol, char *nan.e , 



60 
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int i, j, 
FILE *file? 

file = FOPEN( name, "at" ); 

fort >:«=0; VKnsli; k++ ) 
( 
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for( i-0; J.<nrow; i++ ) 
( 

for( j«0; j<ncol? j++ ) 
{ 

fprintf( file, "%g\t", data (k*nrow*ncol+i*ncol+j3 ); 

} 

**fprintf( file, "\n" ); 

) 

fpL-intf( file, »\n" ); 

1 

FCLOSEl file ); 

) 

IS__BIG_ENDIAN ( } : Testen, ob Big-Endian- oder Little^-Endian-System. 

int IS_BIG_ENDIAN< void ) 
{ 

char *chr; 
short sht; 
int rtn; 

sht = 0x0100; 

chr e (char *)tsht; 

rtn = *chr; 

return rtn; 

) 

SWAP_TWO<) : Vertavischt zwei Bytes <nur intern). 



void SWAP_TWO( char *chrl, char ♦chr 2 ) 
{ 

char tmp; 

tmp «= *chrl; 
*chrl «= *chr2; 
*chr2 « tmp; 

) 

SWAPC): Aendert die Byte-Reihenfolge einer Variablen. 

****-^****«^*************** 1^*4- / 

void SWAP( void *ptr, int size ) 
{ 

char "♦'chr; 
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chr = (char "^)ptr; 

switch { size ) 

case 2 : 

SWAP_TWO( chr, chr+1 ); 
break; 

case 4: 50 
SWAP_TWO{ chr, chr+3 ); 
SWAP_TWO( chr+1, chr+2 ); 
break; 

case &: 

SWAP_TWO( chr, chr+7 ); 55 
SWAP_TWO( chr+1, chr + 6 ); 
SWAP_TWO( chr+2, chr+5 ); 

SWAP_TWO( chr+3, chr +4 ); ^ 
break; 

default: 60 
ERROR( "SWAPO: Wrong size" ); 

) 

) 

[0034] Der grundsaizliche Aufbau der Vorrichtung wird anhand der Figur verdeutiichi. 65 
[0035] Ein Objekt 1 befindet sich naher an einem Multi-Loch-Kollimator 3 als die Delektoroberflache 2. Der Multi- 
Lcx;h-Kollimator 3 weistLocher 4 auf, die von beiden Seiten irichterfbmiig in den KoUimator 3 einmunden, um so einen 
Durchgang von schrag auftreffenden (Tammaquanten durch die Locher zu ermoglichen. Die Spitzen 5 des Multi-Loch- 
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erhaltenen Messergebnis iibereinstimmt. 
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